Graphene, which is a single layer of carbon atoms packed in a honeycomb crystal lattice, has stimulated growing research interests and features wide applications in the fields of electronics and optoelectronics and as energy conversion devices^[1–4] due to its exceptional electronic and optoelectronic properties. Graphene also possesses ultrahigh in-plane thermal conductivity (in the range of 2000–5000 W/mK^[5), which is extremely useful for heat dissipation in next-generation nanodevices. The ultrahigh in-plane thermal conductivity of graphene is attributed to covalent sp² bonding between carbon atoms.^[6 Phonons can transport through pure graphene nanosheets almost without any disturbance due to the large phonon mean free path (∼ 775 nm^[7]). Theoretical and experimental studies showed that the thermal conductivity of graphene shows sensitivity to small perturbations, such as strain,^[8] defects,^[9,10] substrate effect,^[11] and stiffness,^[12] thus allowing the possible modulation of the thermal conductivity of graphene.

The renewed burst in heat transport properties of graphene also presents important implications for other technologies, of which an important example is thermoelectrics. The high thermoelectric conversion efficiency of materials requires simultaneous high electrical conductivity and Seebeck coefficient but low thermal conductivity.^[13] However, graphene possesses naturally poor thermoelectric properties because of its very high thermal conductivity and gapless characteristic, which result in extremely small Seebeck coefficients.^[14] Graphene nanoribbons (GNRs) are strips of graphene with a narrow width that is normally less than 50 nm, indicating notable electrical properties, such as high mobility, high conductance, and small bandgap. Recently, theoretical and experimental studies^[15–17] reported that GNRs show good potential as thermoelectric materials due to their small thermal conductivity without degrading their electronic properties compared with pure graphene.

Previous studies have shown that chemical functionalization serves as an effective routine to modulate the electronic and thermal properties of silicon nanostructures, including silicon nanowires^[18–21] and silicene.^[22] Recent studies of the influence of surface functionalization on the thermal properties of graphene/GNRs have also attracted considerable attention. Zhang et al.^[23] and Sun et al.^[24] performed molecular dynamics (MD) studies on tailoring the thermal transport property of graphene through oxygen functionalization. Chien et al.^[25] investigated the thermal conductivity of GNRs functionalized by chemical attachment of methyl and phenyl groups by using MD simulations and discovered that reduced thermal conductivity strongly depends on the degree of functionalization. Wang et al.^[26] also theoretically reported the effect of surface functionalization on the thermal conductivity of graphene–polymer nanocomposites. In addition, numerous studies^[27–30] have investigated the effect of alkyl functional groups on the structure and surface chemistry properties of graphene/GNRs, but the influence of alkyl chain functionalization on the thermal conductivity of GNRs has not been well understood so far.

In this study, we performed reverse non-equilibrium MD (RNEMD) simulation to investigate the thermal conductivity of alkyl group-functionalized zigzag GNRs (ZGNRs). The effect of alkyl group functionalization, i.e., grafting hydrocarbon chains to ZGNRs with covalent bonds, on phonon thermal conductivity was investigated and compared with pristine ZGNRs. The effects of length and concentration of alkyl chains on the thermal conductivity of ZGNRs were considered, and the underlying mechanism was explained based on phonon spectrum analysis. Our study reveals the significant influence of alkyl chain functionalization on phonon thermal transport in ZGNRs, benefiting the thermoelectric applications of GNR-based materials.

Phonon thermal conductivity can be predicted by either equilibrium or non-equilibrium MD (NEMD) methods. Both methods feature their own advantages and disadvantages. Equilibrium MD calculates thermal conductivity using Einstein^[31] or Green–Kubo relations.^[13] In this approach, the system will remain at equilibrium conditions. By contrast, in the NEMD method, a temperature gradient is applied in the system, and the resulting heat flux is measured. The simulation must run long enough to reach the non-equilibrium steady state, where constant heat flux and temperature profile are observed. Alternatively, a heat flux can also be directly imposed on the system by adding or subtracting kinetic energies, which determine the temperature gradient induced, a condition called RNEMD. Based on the framework of RNEMD, Müller-Plathe^[32] proposed a simple and fast converging scheme, which offers both conservation of total energy and conservation of total linear momentum.

Here, the Müller-Plathe-based RNEMD method was performed to calculate the thermal conductivity of ZGNR–alkyl chain systems using LAMMPS simulation software.^[33] First, a simple ZGNR sheet was sketched. Then, the alkyl chains with different lengths and functionalization concentrations were placed on a ZGNR basal plane. In MD simulations, the periodic boundary condition was applied along the length and width of the ZGNR sheet (x and y directions, respectively). A time step of 0.5 fs was considered for integration of the atomic equations of motion. Carbon–carbon and carbon–hydrogen interactions were described using the adaptive intermolecular reactive empirical bond order (AIREBO) potential, which is an extension of the second-generation reactive bond order (REBO) potential,^[34] and has been successfully implemented to study the thermal and mechanical properties of carbon-based systems, including graphene.^[25,35–37]

The initial configuration was relaxed by MD simulation in the isothermal-isobaric (NPT) ensemble at room temperature for 2 × 10⁶ time steps (1 ns). Then the structure was first equilibrated using the canonical (NVT) ensemble simulation for 1 ns under the Nose–Hoover thermostat^[38] with a temperature of 300 K. The average length of carbon–carbon bond in ZGNR and that in alkyl chains remained at around 1.42 Å and 1.54 Å, respectively, which are very close to the first-principle calculations.^[39,40] After subsequent 1 ns simulation with the microcanonical (NVE) ensemble, the system keeps the total energy constant. Before RNEMD calculations, all the GNR structures were fully optimized. Next, the Müller-Plathe MD simulation was carried out in the NVE ensemble to calculate phonon thermal conductivity. The simulation system was divided into N slabs along the heat flux direction (x direction). Two slabs in the left and right edges were considered cold slabs, and the middle slab was considered a hot slab, as shown in Fig. 1(a). The thickness of the hot region was twice that of the cold region to ensure that the number of atoms in the hot region was the same as that in the cold region. Heat flux J_x was created by exchanging the momentum between the hottest atom velocity in the cold (sink) region and the coldest atom velocity in the hot (source) region every swap time interval of 50 fs, inducing a temperature gradient in the system of interest. Under the NVE ensemble, we ran 2 × 10⁶ time steps (1 ns) to allow the system to reach a steady state. After N_swap exchange was performed, heat flux J_x was given by the sum of exchanged kinetic energy per unit time and area as

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (a) Schematic of the simulation model for the Müller–Plathe RNEMD method. Heat flux was established between the hot and cold regions. Periodic boundary conditions were applied in the x and y directions. (b) Typical temperature profile in a ZGNR along the x direction at 300 K and the temperature gradient calculated by the linear response region of the temperature profile.

We also analyzed the phonon density of states (PDOS) along the direction of the temperature gradient to gain insights into the mechanism of the thermal behavior of GNRs. PDOS as a function of the phonon angular frequency ω can be obtained from fast Fourier transform of the velocity autocorrelation function^[42] as follows:

The accuracy of the used simulation methods was validated by considering a pristine ZGNR with 16 nm length and 4 nm width as a model system and calculating the phonon thermal conductivity of this model system, which is divided into two sections in the RNEMD method. By assuming the thickness of the ZGNR to be 0.142 nm,^[43] the calculated phonon thermal conductivity of the model system measured 119.02 ± 3.91 W/mK at 300 K, showing good agreement with previously AIREBO potential-based simulated values of 135.01 ± 13.38 W/mK for a 17 nm-long ZGNR^[44] and 140 W/mK for a 20 nm-long ZGNR.^[45] In our test practice, we observed the considerable sensitivity of thermal conductivity of the studied ZGNR to ZGNR length due to acoustic phonon scattering. For example, when the length of ZGNR was doubled (i.e., 32 nm), the calculated thermal conductivity increased to about twice that of the original short ZGNR (i.e., 220.57 ± 5.68 W/mK). The experimentally reported thermal conductivity values of micrometer-sized graphene ranged from 2000 W/mK to 5000 W/mK,^[5,7] which are significantly larger than our MD simulation results. This difference can be attributed to the length dependence of thermal transport in nanosheets.^[18,43,44] The mean free path of phonons in graphene approximated 775 nm,^[7] whereas the length of our simulation cell approximated 16 nm. Therefore, phonon-boundary scattering dominated and caused remarkable reduction of the thermal conductivity in GNRs.^[36]

In the present work, we limited our attention to the effect of functional groups rather than the absolute value of phonon thermal conductivity. Therefore, we used relative thermal conductivity, which is the ratio between the thermal conductivities of functionalized and pristine GNRs to avoid the effects of different thickness assumptions and size length and boundary.

Figure 2(a) shows the phonon thermal conductivity of GNR monofunctionalized by an alkyl chain with different lengths. The phonon thermal conductivity of ZGNR decreased rapidly (by about 18%) after a methyl group (n = 1) was integrated. With the further increase in the number of carbon atoms in the alkyl chain, the phonon thermal conductivity of ZGNR decreased gradually. For example, pentadecyl chain (n = 15) monofunctionalization lowered the phonon thermal conductivity of the pristine GNR without surface functionalization by about 25%. Figure 2(b) shows the influence of functionalization concentration on the phonon thermal conductivity of methyl group-functionalized GNR. The concentration of functionalization (ρ) is defined as the ratio of the number of sp³-hybridized carbon atoms to the number of total carbon atoms in GNR. In the present study, the concentration of functionalization (ρ) varied from 0% to 0.1% for methyl group random functionalization. We observed that the phonon thermal conductivity of GNR was much more sensitive to alkyl chain concentration than alkyl chain length. The simulations predicted a decrease in phonon thermal conductivity by as much as 18% as one carbon atom in the GNR was functionalized. Functionalization of 0.01% to 0.1% continually degraded the phonon thermal conductivity from ∼80% to ∼50% compared with the studied pristine GNR. Particularly, when the concentration of functionalization was increased to above 0.05%, the thermal conductivity became less sensitive to the functional group coverage. Similar phenomena were also observed in phenyl-functionalized GNRs.^[25]

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. The calculated phonon thermal conductivity of the functionalized ZGNR as functions of (a) alkyl carbon atom number (n) for monofunctionalization and (b) functionalization concentration (ρ) of methyl group, compared with the pristine ZGNR. Here, the studied ZGNR measured 16 nm in length and 4 nm in width.

Phonon scattering was investigated to understand the role of alkyl groups on phonon thermal transport by calculating and comparing the x direction PDOS of the different segments at the two ends of the nanoribbons (Fig. 3). The mismatch between two spectra indicates that phonon vibrations have scattered along the nanosheets.^[46] For pristine ZGNR, our calculation showed the location of the highest frequency carbon–carbon optical mode (corresponding to the G peak) neared 50 THz (see Fig. 3(a)), which is in good agreement with previously AIREBO potential-based simulated results^[36,45] and experimental result (47 THz^[47]). The PDOS of the left and right segments were almost identical throughout the entire frequency range, indicating that phonons easily traveled through pristine GNR in the x direction with no disturbance. Our study also shows that functionalization of GNR with alkyl functional groups led to a large mismatch between the spectra of the left and right segments. For example, as shown in Figs. 3(b)–3(d), with increasing methyl group concentration from 0.025% to 0.05% and to 0.1%, the G peak of the right segment was effectively suppressed because the corresponding intensity considerably decreased compared with that of the left segment. However, the low-frequency acoustic modes between two segments caused no significant change under functionalization of methyl groups.

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. Calculated PDOS along the x direction of the atoms in the left and right segments for pristine ZGNR sheet and methyl group-functionalized ZGNR with different functionalization concentrations (ρ), respectively: (a) ρ = 0, (b) ρ = 0.025%, (c) ρ = 0.05%, (d) ρ = 0.1%. The top panel shows the left (red color) and right segments (blue color) at the two edges of the nanosheet. The mismatch between the two spectra indicates that phonon vibrations have scattered along the nanosheets.

When the alkyl chain is connected to the ZGNR sheet, the atoms on which the alkyl chains are covalently bonded contain hybridizations that change from sp² to sp³. The introduction of sp³ carbons creates structural defect sites in lattice for the phonons to interact with and scatter off,^[48] thus reducing the mean free path of the phonons. The functionalization-induced phonon spectra mismatch of the two segments indicates that phonons cannot travel easily through the nanosheet owing to remarkable phonon scattering. Therefore, increasing the length and/or concentration of alkyl groups can induce the enhancement of phonon–structural defect scattering in ZGNR, thus resulting in large decrement of phonon thermal conductivity of the ZGNR. In addition, the perturbation of alkyl chain in the dynamics can apply local stress field to the ZGNR sheet and thus induce phonon localization in graphene nanosheets due to local strain,^[49] which can further decrease the phonon thermal conductivity.

We have investigated the influence of alkyl group functionalization on phonon thermal conductivity and phonon spectrum of ZGNR by using the RNEMD method. Our study showed the reduction of phonon thermal conductivity in ZGNR induced by functionalization of alkyl groups. The concentration of alkyl group can be a major parameter for modulating the thermal conductivity of ZGNR due to phonon scattering. The length of alkyl chain can also influence the phonon thermal conductivity of ZGNR through phonon localization. Phonon spectra analysis showed that functionalization of ZGNR with alkyl functional groups led to a large mismatch between the spectra of the left and right segments. The spectra mismatch increased with increasing chain length and concentration of the alkyl functional groups, resulting in decreased phonon thermal conductivity. Our study can provide a useful reference for designing graphene–polymer nanocomposites with precisely controlled thermal property for thermal management and thermoelectric applications.